Figure3D

Load Data

library(gdata)
dataset <- read.delim("raw_data/Figure3D.txt", stringsAsFactors = FALSE)

dataset$genotype <- factor(gsub(" ","",dataset$genotype))
dataset$experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each= length(levels(dataset$genotype))))

dataset$UID <- factor(paste(dataset$experiment, dataset$genotype))

kable(dataset, row.names = F)

genotype	NT	X2h	X4h	X6h	experiment	UID
WT	8.366667	35.60000	29.56667	21.23333	exp1	exp1 WT
PARP1KO	6.700000	57.23333	46.26667	45.96667	exp1	exp1 PARP1KO
ALC1KO	6.533333	56.96667	59.76667	55.66667	exp1	exp1 ALC1KO
ALC1KOPARP1KO	6.400000	55.26667	43.00000	46.40000	exp1	exp1 ALC1KOPARP1KO
WT	7.300000	36.13333	29.73333	23.00000	exp2	exp2 WT
PARP1KO	7.966667	48.96667	49.60000	43.93333	exp2	exp2 PARP1KO
ALC1KO	10.600000	54.16667	59.83333	57.93333	exp2	exp2 ALC1KO
ALC1KOPARP1KO	8.300000	53.43333	46.23333	42.43333	exp2	exp2 ALC1KOPARP1KO
WT	7.233333	34.06667	35.43333	20.66667	exp3	exp3 WT
PARP1KO	7.233333	39.76667	41.26667	33.56667	exp3	exp3 PARP1KO
ALC1KO	7.233333	54.46667	52.10000	55.13333	exp3	exp3 ALC1KO
ALC1KOPARP1KO	6.433333	42.43333	42.90000	35.96667	exp3	exp3 ALC1KOPARP1KO

library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "treatment", value.name = "Outcome")

dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT")

dataset$time <- gsub("NT","0",dataset$treatment)
dataset$time <- gsub("X|h","",dataset$time)
dataset$time <- as.integer(dataset$time)



kable(dataset, row.names = F)

genotype	experiment	UID	treatment	Outcome	time
WT	exp1	exp1 WT	NT	8.366667	0
PARP1KO	exp1	exp1 PARP1KO	NT	6.700000	0
ALC1KO	exp1	exp1 ALC1KO	NT	6.533333	0
ALC1KOPARP1KO	exp1	exp1 ALC1KOPARP1KO	NT	6.400000	0
WT	exp2	exp2 WT	NT	7.300000	0
PARP1KO	exp2	exp2 PARP1KO	NT	7.966667	0
ALC1KO	exp2	exp2 ALC1KO	NT	10.600000	0
ALC1KOPARP1KO	exp2	exp2 ALC1KOPARP1KO	NT	8.300000	0
WT	exp3	exp3 WT	NT	7.233333	0
PARP1KO	exp3	exp3 PARP1KO	NT	7.233333	0
ALC1KO	exp3	exp3 ALC1KO	NT	7.233333	0
ALC1KOPARP1KO	exp3	exp3 ALC1KOPARP1KO	NT	6.433333	0
WT	exp1	exp1 WT	X2h	35.600000	2
PARP1KO	exp1	exp1 PARP1KO	X2h	57.233333	2
ALC1KO	exp1	exp1 ALC1KO	X2h	56.966667	2
ALC1KOPARP1KO	exp1	exp1 ALC1KOPARP1KO	X2h	55.266667	2
WT	exp2	exp2 WT	X2h	36.133333	2
PARP1KO	exp2	exp2 PARP1KO	X2h	48.966667	2
ALC1KO	exp2	exp2 ALC1KO	X2h	54.166667	2
ALC1KOPARP1KO	exp2	exp2 ALC1KOPARP1KO	X2h	53.433333	2
WT	exp3	exp3 WT	X2h	34.066667	2
PARP1KO	exp3	exp3 PARP1KO	X2h	39.766667	2
ALC1KO	exp3	exp3 ALC1KO	X2h	54.466667	2
ALC1KOPARP1KO	exp3	exp3 ALC1KOPARP1KO	X2h	42.433333	2
WT	exp1	exp1 WT	X4h	29.566667	4
PARP1KO	exp1	exp1 PARP1KO	X4h	46.266667	4
ALC1KO	exp1	exp1 ALC1KO	X4h	59.766667	4
ALC1KOPARP1KO	exp1	exp1 ALC1KOPARP1KO	X4h	43.000000	4
WT	exp2	exp2 WT	X4h	29.733333	4
PARP1KO	exp2	exp2 PARP1KO	X4h	49.600000	4
ALC1KO	exp2	exp2 ALC1KO	X4h	59.833333	4
ALC1KOPARP1KO	exp2	exp2 ALC1KOPARP1KO	X4h	46.233333	4
WT	exp3	exp3 WT	X4h	35.433333	4
PARP1KO	exp3	exp3 PARP1KO	X4h	41.266667	4
ALC1KO	exp3	exp3 ALC1KO	X4h	52.100000	4
ALC1KOPARP1KO	exp3	exp3 ALC1KOPARP1KO	X4h	42.900000	4
WT	exp1	exp1 WT	X6h	21.233333	6
PARP1KO	exp1	exp1 PARP1KO	X6h	45.966667	6
ALC1KO	exp1	exp1 ALC1KO	X6h	55.666667	6
ALC1KOPARP1KO	exp1	exp1 ALC1KOPARP1KO	X6h	46.400000	6
WT	exp2	exp2 WT	X6h	23.000000	6
PARP1KO	exp2	exp2 PARP1KO	X6h	43.933333	6
ALC1KO	exp2	exp2 ALC1KO	X6h	57.933333	6
ALC1KOPARP1KO	exp2	exp2 ALC1KOPARP1KO	X6h	42.433333	6
WT	exp3	exp3 WT	X6h	20.666667	6
PARP1KO	exp3	exp3 PARP1KO	X6h	33.566667	6
ALC1KO	exp3	exp3 ALC1KO	X6h	55.133333	6
ALC1KOPARP1KO	exp3	exp3 ALC1KOPARP1KO	X6h	35.966667	6

Plot Data

library(ggplot2)

# quadratic
ggplot(dataset, aes(x=time, y=Outcome)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=genotype)) +
        geom_point(aes(colour=genotype, shape=experiment), size=2) +        
        facet_grid(. ~ genotype) +
        xlab(label = "Time post 2 Gy irradiation") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#000080","#808080","#800000"))

# cubic
ggplot(dataset, aes(x=time, y=Outcome)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE, aes(colour=genotype)) +
        geom_point(aes(colour=genotype, shape=experiment), size=2) +        
        facet_grid(. ~ genotype) +
        xlab(label = "Time post 2 Gy irradiation") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#000080","#808080","#800000"))

dataplot <- aggregate(Outcome ~ genotype+treatment, data = dataset, mean)
dataplot$se <- aggregate(Outcome ~ genotype+treatment, data = dataset, FUN = function(x){ sd(x)/sqrt(length(x))})$Outcome
dataplot$treatment <- factor(gsub("X|.nM","", dataplot$treatment)) 

dataplot$genotype <- factor(dataplot$genotype, levels = c("WT","PARP1KO","ALC1KO","ALC1KOPARP1KO"))
dataplot$treatment <- relevel(dataplot$treatment, ref = "NT")


ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_bar(stat="identity", position=position_dodge()) +
        geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9)) +
        ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
        xlab(label = "Time post 2 Gy irradiation") +
        ylim(0, 70) +
        scale_fill_manual(values=c("#000000","#800000","#000080","#808080"))

library(Cairo)

cairo_pdf("Figure3D_v1.pdf", width = 7, height = 4, family = "Arial")

ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) + 
        theme_bw() +
        theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
              axis.line = element_line(colour = "black"), text = element_text(size=14),
              panel.border = element_blank(), panel.background = element_blank()) +
        geom_bar(stat="identity", position=position_dodge()) +
        geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9)) +
        ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
        xlab(label = "Time post 2 Gy irradiation (h)") +
        ylim(0, 70) +
        scale_fill_manual(values=c("#000000","#800000","#000080","#808080"))

dev.off()

## quartz_off_screen 
##                 2

dataset$genotype <- factor(dataset$genotype, levels = c("WT","PARP1KO","ALC1KO","ALC1KOPARP1KO"))
dataset$treatment <- factor(gsub("X","", dataset$treatment)) 
dataset$treatment <- relevel(dataset$treatment, ref = "NT")



cairo_pdf("Figure3D_v2.pdf", width = 7, height = 4, family = "Arial")

ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) + 
        theme_bw() +
        theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
              axis.line = element_line(colour = "black"), text = element_text(size=14),
              panel.border = element_blank(), panel.background = element_blank()) +
        geom_bar(stat="identity", position=position_dodge(), colour="black") +
        ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
        xlab(label = "Time post 2 Gy irradiation (h)") +
        ylim(0, 70) +
        scale_fill_manual(values=c("#00000080","#80000080","#00008080","#80808080")) +
        geom_jitter(data = dataset, cex=1, position=position_dodge(0.9),  aes(colour=genotype)) +
        scale_color_manual(values=c("#000000","#800000","#000080","#808080")) +
        geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9)) 


cairo_pdf("Figure3D_v3.pdf", width = 7, height = 4, family = "Arial")

ggplot(dataset, aes(x=treatment, y=Outcome, color=genotype)) +
        theme_bw() +
        theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
              axis.line = element_line(colour = "black"), text = element_text(size=14),
              panel.border = element_blank(), panel.background = element_blank()) +
        geom_jitter(cex=1, position=position_dodge(1)) +
        stat_summary(fun.data=mean_se, fun.args = list(mult=1), geom="errorbar", width=0.25, aes(colour=genotype), position=position_dodge(1)) +
        stat_summary(fun.y=mean, geom="crossbar", width=0.5, aes(colour=genotype), position=position_dodge(1)) +
        ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
        xlab(label = "Time post 2 Gy irradiation (h)") +
        ylim(0, 70) +
        scale_fill_manual(values=c("#00000080","#80000080","#00008080","#80808080")) +
        scale_color_manual(values=c("#000000","#800000","#000080","#808080"))


dev.off()

## quartz_off_screen 
##                 2

Linear Model

library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)

fit1 <- lm(Outcome ~ experiment+time*genotype, data = dataset)
print(summary(fit1))

## 
## Call:
## lm(formula = Outcome ~ experiment + time * genotype, data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.090 -11.403  -2.183  11.781  25.076 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                19.66444    7.29796   2.695   0.0104 *
## experimentexp2             -0.08542    4.92028  -0.017   0.9862  
## experimentexp3             -4.06458    4.92028  -0.826   0.4139  
## time                        1.91556    1.79663   1.066   0.2931  
## genotypePARP1KO             2.63111    9.50689   0.277   0.7835  
## genotypeALC1KO              3.95889    9.50689   0.416   0.6794  
## genotypeALC1KOPARP1KO       2.88556    9.50689   0.304   0.7631  
## time:genotypePARP1KO        3.01556    2.54082   1.187   0.2427  
## time:genotypeALC1KO         5.40444    2.54082   2.127   0.0400 *
## time:genotypeALC1KOPARP1KO  2.95111    2.54082   1.161   0.2527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.92 on 38 degrees of freedom
## Multiple R-squared:  0.5491, Adjusted R-squared:  0.4423 
## F-statistic: 5.141 on 9 and 38 DF,  p-value: 0.0001444

cat("AIC: ", AIC(fit1))

## AIC:  399.7809

simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(Outcome ~ experiment+poly(time,2)*genotype, data = dataset)
print(summary(fit2))

## 
## Call:
## lm(formula = Outcome ~ experiment + poly(time, 2) * genotype, 
##     data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4611  -4.2647  -0.0672   3.4001  13.5978 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                           25.41111    2.28095  11.141 6.84e-13 ***
## experimentexp2                        -0.08542    2.28095  -0.037   0.9703    
## experimentexp3                        -4.06458    2.28095  -1.782   0.0837 .  
## poly(time, 2)1                        29.67566   12.90300   2.300   0.0277 *  
## poly(time, 2)2                       -65.08662   12.90300  -5.044 1.51e-05 ***
## genotypePARP1KO                       11.67778    2.63381   4.434 9.19e-05 ***
## genotypeALC1KO                        20.17222    2.63381   7.659 6.66e-09 ***
## genotypeALC1KOPARP1KO                 11.73889    2.63381   4.457 8.59e-05 ***
## poly(time, 2)1:genotypePARP1KO        46.71679   18.24760   2.560   0.0151 *  
## poly(time, 2)2:genotypePARP1KO       -14.43376   18.24760  -0.791   0.4344    
## poly(time, 2)1:genotypeALC1KO         83.72529   18.24760   4.588 5.83e-05 ***
## poly(time, 2)2:genotypeALC1KO        -18.16729   18.24760  -0.996   0.3265    
## poly(time, 2)1:genotypeALC1KOPARP1KO  45.71842   18.24760   2.505   0.0172 *  
## poly(time, 2)2:genotypeALC1KOPARP1KO -14.20282   18.24760  -0.778   0.4418    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.452 on 34 degrees of freedom
## Multiple R-squared:  0.9133, Adjusted R-squared:  0.8801 
## F-statistic: 27.55 on 13 and 34 DF,  p-value: 2.899e-14

cat("AIC: ", AIC(fit2))

## AIC:  328.6398

simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)
print(summary(fit3))

## 
## Call:
## lm(formula = Outcome ~ experiment + poly(time, 3) * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2076 -1.9667  0.0132  1.7253  7.1944 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                           25.41111    1.21180  20.970  < 2e-16 ***
## experimentexp2                        -0.08542    1.21180  -0.070 0.944274    
## experimentexp3                        -4.06458    1.21180  -3.354 0.002170 ** 
## poly(time, 3)1                        29.67566    6.85499   4.329 0.000154 ***
## poly(time, 3)2                       -65.08662    6.85499  -9.495 1.51e-10 ***
## poly(time, 3)3                        19.41656    6.85499   2.832 0.008176 ** 
## genotypePARP1KO                       11.67778    1.39927   8.346 2.58e-09 ***
## genotypeALC1KO                        20.17222    1.39927  14.416 5.01e-15 ***
## genotypeALC1KOPARP1KO                 11.73889    1.39927   8.389 2.31e-09 ***
## poly(time, 3)1:genotypePARP1KO        46.71679    9.69443   4.819 3.88e-05 ***
## poly(time, 3)2:genotypePARP1KO       -14.43376    9.69443  -1.489 0.146959    
## poly(time, 3)3:genotypePARP1KO        13.65011    9.69443   1.408 0.169399    
## poly(time, 3)1:genotypeALC1KO         83.72529    9.69443   8.636 1.24e-09 ***
## poly(time, 3)2:genotypeALC1KO        -18.16729    9.69443  -1.874 0.070701 .  
## poly(time, 3)3:genotypeALC1KO         13.13372    9.69443   1.355 0.185606    
## poly(time, 3)1:genotypeALC1KOPARP1KO  45.71842    9.69443   4.716 5.19e-05 ***
## poly(time, 3)2:genotypeALC1KOPARP1KO -14.20282    9.69443  -1.465 0.153310    
## poly(time, 3)3:genotypeALC1KOPARP1KO  22.06740    9.69443   2.276 0.030131 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.427 on 30 degrees of freedom
## Multiple R-squared:  0.9784, Adjusted R-squared:  0.9662 
## F-statistic: 79.96 on 17 and 30 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit3))

## AIC:  269.9136

simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

Compare Results

ICtab(fit1,fit2,fit3,
      base=T)

##      AIC   dAIC  df
## fit3 269.9   0.0 19
## fit2 328.6  58.7 15
## fit1 399.8 129.9 11

Final Result

fit <- fit3

output <- coef(summary(fit))
output <- output[grep("time", rownames(output)),]

rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )


# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
time1 in WT	29.67566	6.854994	4.329057	0.0001535
time2 in WT	-65.08662	6.854994	-9.494774	0.0000000
time3 in WT	19.41656	6.854994	2.832469	0.0081758
time1: WT vs. PARP1KO	46.71679	9.694426	4.818933	0.0000388
time2: WT vs. PARP1KO	-14.43376	9.694426	-1.488872	0.1469593
time3: WT vs. PARP1KO	13.65011	9.694426	1.408038	0.1693989
time1: WT vs. ALC1KO	83.72529	9.694426	8.636437	0.0000000
time2: WT vs. ALC1KO	-18.16729	9.694426	-1.873993	0.0707005
time3: WT vs. ALC1KO	13.13372	9.694426	1.354770	0.1856062
time1: WT vs. ALC1KOPARP1KO	45.71842	9.694426	4.715949	0.0000519
time2: WT vs. ALC1KOPARP1KO	-14.20282	9.694426	-1.465050	0.1533098
time3: WT vs. ALC1KOPARP1KO	22.06740	9.694426	2.276298	0.0301313

write.table(output, file = "Figure3D_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1KO")

fit <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)

output <- coef(summary(fit))
output <- output[grep("time", rownames(output)),]

rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )


# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
time1 in ALC1KO	113.4009524	6.854994	16.5428227	0.0000000
time2 in ALC1KO	-83.2539088	6.854994	-12.1450008	0.0000000
time3 in ALC1KO	32.5502734	6.854994	4.7484028	0.0000474
time1: ALC1KO vs. WT	-83.7252933	9.694426	-8.6364367	0.0000000
time2: ALC1KO vs. WT	18.1672885	9.694426	1.8739933	0.0707005
time3: ALC1KO vs. WT	-13.1337169	9.694426	-1.3547700	0.1856062
time1: ALC1KO vs. PARP1KO	-37.0085075	9.694426	-3.8175039	0.0006289
time2: ALC1KO vs. PARP1KO	3.7335317	9.694426	0.3851215	0.7028645
time3: ALC1KO vs. PARP1KO	0.5163978	9.694426	0.0532675	0.9578719
time1: ALC1KO vs. ALC1KOPARP1KO	-38.0068766	9.694426	-3.9204877	0.0004749
time2: ALC1KO vs. ALC1KOPARP1KO	3.9644719	9.694426	0.4089434	0.6854860
time3: ALC1KO vs. ALC1KOPARP1KO	8.9336816	9.694426	0.9215277	0.3641269

write.table(output, file = "Figure3D_Stats_Ref_ALC1KO.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

ANOVA

fit3a <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)
fit3b <- lm(Outcome ~ experiment+poly(time,3)+genotype, data = dataset)

anova(fit3a, fit3b)

## Analysis of Variance Table
## 
## Model 1: Outcome ~ experiment + poly(time, 3) * genotype
## Model 2: Outcome ~ experiment + poly(time, 3) + genotype
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     30  352.43                                  
## 2     39 1343.95 -9   -991.51 9.3778 1.297e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

New analysis

• fit model for each time point independently

dataset$genotype <- relevel(dataset$genotype, ref = "WT")


for(i in 2:length(levels(dataset$treatment))){
        
        dataset_sub <- dataset[dataset$treatment == levels(dataset$treatment)[i],]
        
        for(g in seq_along(levels(dataset_sub$genotype))){
                
                
                cat("Time: ", levels(dataset$treatment)[i], "\n",
                    "Genotype: ", levels(dataset$genotype)[g], "\n")
                
                
                dataset_sub$genotype <- relevel(dataset_sub$genotype, ref = levels(dataset_sub$genotype)[g])

                fit_sub <- lm(Outcome ~ genotype, data = dataset_sub)
                
                print(summary(fit_sub))
                cat("AIC: ", AIC(fit_sub))
                
                
                simres <- simulateResiduals(fittedModel = fit_sub)
                plot(simres)
                
                if(i == 2 & g == 1){
                        output <- coef(summary(fit_sub))
                        output <- output[grep("genotype", rownames(output)),]
                        rownames(output) <- paste0(gsub("genotype",  paste0(" ",levels(dataset$genotype)[g], " vs. "), rownames(output)),
                                                   " in Time ", levels(dataset$treatment)[i])
                } else {
                        outtmp <- coef(summary(fit_sub))
                        outtmp <- outtmp[grep("genotype", rownames(outtmp)),]
                        rownames(outtmp) <- paste0(gsub("genotype",  paste0(" ",levels(dataset$genotype)[g], " vs. "), rownames(outtmp)),
                                                   " in Time ", levels(dataset$treatment)[i])
                        
                        output <- rbind(output,outtmp)
                        
                        if(g == length(levels(dataset_sub$genotype)) & i < length(levels(dataset$treatment))){
                                
                                output <- rbind(output, " ", colnames(output))
                        }
                }
                
        }
        
        
}

## Time:  2h 
##  Genotype:  WT 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8889 -1.0750  0.3222  2.0889  8.5778 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             35.267      3.266  10.797 4.78e-06 ***
## genotypeALC1KO          19.933      4.619   4.315  0.00256 ** 
## genotypePARP1KO         13.389      4.619   2.898  0.01994 *  
## genotypeALC1KOPARP1KO   15.111      4.619   3.271  0.01134 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.6137 
## F-statistic: 6.826 on 3 and 8 DF,  p-value: 0.01349
## 
## AIC:  80.78097

## Time:  2h 
##  Genotype:  ALC1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8889 -1.0750  0.3222  2.0889  8.5778 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             55.200      3.266  16.899 1.52e-07 ***
## genotypeWT             -19.933      4.619  -4.315  0.00256 ** 
## genotypePARP1KO         -6.544      4.619  -1.417  0.19431    
## genotypeALC1KOPARP1KO   -4.822      4.619  -1.044  0.32705    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.6137 
## F-statistic: 6.826 on 3 and 8 DF,  p-value: 0.01349
## 
## AIC:  80.78097

## Time:  2h 
##  Genotype:  PARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8889 -1.0750  0.3222  2.0889  8.5778 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             48.656      3.266  14.896 4.07e-07 ***
## genotypeALC1KO           6.544      4.619   1.417   0.1943    
## genotypeWT             -13.389      4.619  -2.898   0.0199 *  
## genotypeALC1KOPARP1KO    1.722      4.619   0.373   0.7190    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.6137 
## F-statistic: 6.826 on 3 and 8 DF,  p-value: 0.01349
## 
## AIC:  80.78097

## Time:  2h 
##  Genotype:  ALC1KOPARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8889 -1.0750  0.3222  2.0889  8.5778 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       50.378      3.266  15.423 3.11e-07 ***
## genotypePARP1KO   -1.722      4.619  -0.373   0.7190    
## genotypeALC1KO     4.822      4.619   1.044   0.3270    
## genotypeWT       -15.111      4.619  -3.271   0.0113 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.6137 
## F-statistic: 6.826 on 3 and 8 DF,  p-value: 0.01349
## 
## AIC:  80.78097

## Time:  4h 
##  Genotype:  WT 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1333 -1.8861 -0.2444  2.5500  3.8889 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             31.578      2.084  15.154 3.56e-07 ***
## genotypeALC1KO          25.656      2.947   8.706 2.36e-05 ***
## genotypePARP1KO         14.133      2.947   4.796  0.00136 ** 
## genotypeALC1KOPARP1KO   12.467      2.947   4.230  0.00288 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared:  0.9049, Adjusted R-squared:  0.8693 
## F-statistic: 25.39 on 3 and 8 DF,  p-value: 0.0001931
## 
## AIC:  69.99344

## Time:  4h 
##  Genotype:  ALC1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1333 -1.8861 -0.2444  2.5500  3.8889 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             57.233      2.084  27.465 3.33e-09 ***
## genotypeWT             -25.656      2.947  -8.706 2.36e-05 ***
## genotypePARP1KO        -11.522      2.947  -3.910  0.00448 ** 
## genotypeALC1KOPARP1KO  -13.189      2.947  -4.475  0.00207 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared:  0.9049, Adjusted R-squared:  0.8693 
## F-statistic: 25.39 on 3 and 8 DF,  p-value: 0.0001931
## 
## AIC:  69.99344

## Time:  4h 
##  Genotype:  PARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1333 -1.8861 -0.2444  2.5500  3.8889 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             45.711      2.084  21.936 1.97e-08 ***
## genotypeALC1KO          11.522      2.947   3.910  0.00448 ** 
## genotypeWT             -14.133      2.947  -4.796  0.00136 ** 
## genotypeALC1KOPARP1KO   -1.667      2.947  -0.566  0.58721    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared:  0.9049, Adjusted R-squared:  0.8693 
## F-statistic: 25.39 on 3 and 8 DF,  p-value: 0.0001931
## 
## AIC:  69.99344

## Time:  4h 
##  Genotype:  ALC1KOPARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1333 -1.8861 -0.2444  2.5500  3.8889 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       44.044      2.084  21.136 2.64e-08 ***
## genotypePARP1KO    1.667      2.947   0.566  0.58721    
## genotypeALC1KO    13.189      2.947   4.475  0.00207 ** 
## genotypeWT       -12.467      2.947  -4.230  0.00288 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared:  0.9049, Adjusted R-squared:  0.8693 
## F-statistic: 25.39 on 3 and 8 DF,  p-value: 0.0001931
## 
## AIC:  69.99344

## Time:  6h 
##  Genotype:  WT 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5889 -1.0028  0.2167  1.9611  4.8111 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             21.633      2.511   8.616 2.55e-05 ***
## genotypeALC1KO          34.611      3.551   9.747 1.03e-05 ***
## genotypePARP1KO         19.522      3.551   5.498 0.000575 ***
## genotypeALC1KOPARP1KO   19.967      3.551   5.623 0.000497 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared:  0.9231, Adjusted R-squared:  0.8942 
## F-statistic: 31.99 on 3 and 8 DF,  p-value: 8.358e-05
## 
## AIC:  74.46729

## Time:  6h 
##  Genotype:  ALC1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5889 -1.0028  0.2167  1.9611  4.8111 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             56.244      2.511  22.400 1.67e-08 ***
## genotypeWT             -34.611      3.551  -9.747 1.03e-05 ***
## genotypePARP1KO        -15.089      3.551  -4.249  0.00280 ** 
## genotypeALC1KOPARP1KO  -14.644      3.551  -4.124  0.00333 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared:  0.9231, Adjusted R-squared:  0.8942 
## F-statistic: 31.99 on 3 and 8 DF,  p-value: 8.358e-05
## 
## AIC:  74.46729

## Time:  6h 
##  Genotype:  PARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5889 -1.0028  0.2167  1.9611  4.8111 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            41.1556     2.5109  16.391 1.93e-07 ***
## genotypeALC1KO         15.0889     3.5509   4.249 0.002802 ** 
## genotypeWT            -19.5222     3.5509  -5.498 0.000575 ***
## genotypeALC1KOPARP1KO   0.4444     3.5509   0.125 0.903482    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared:  0.9231, Adjusted R-squared:  0.8942 
## F-statistic: 31.99 on 3 and 8 DF,  p-value: 8.358e-05
## 
## AIC:  74.46729

## Time:  6h 
##  Genotype:  ALC1KOPARP1KO 
## 
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5889 -1.0028  0.2167  1.9611  4.8111 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      41.6000     2.5109  16.568 1.78e-07 ***
## genotypePARP1KO  -0.4444     3.5509  -0.125 0.903482    
## genotypeALC1KO   14.6444     3.5509   4.124 0.003325 ** 
## genotypeWT      -19.9667     3.5509  -5.623 0.000497 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared:  0.9231, Adjusted R-squared:  0.8942 
## F-statistic: 31.99 on 3 and 8 DF,  p-value: 8.358e-05
## 
## AIC:  74.46729

# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
WT vs. ALC1KO in Time 2h	19.9333333366666	4.61941688316092	4.31511895133979	0.00256288937859708
WT vs. PARP1KO in Time 2h	13.38888889	4.61941688316092	2.89839372125219	0.0199407600005576
WT vs. ALC1KOPARP1KO in Time 2h	15.11111111	4.61941688316092	3.27121614961496	0.0113355916961444
ALC1KO vs. WT in Time 2h	-19.9333333366666	4.61941688316091	-4.31511895133979	0.00256288937859707
ALC1KO vs. PARP1KO in Time 2h	-6.54444444666667	4.61941688316091	-1.41672523008759	0.194305008686908
ALC1KO vs. ALC1KOPARP1KO in Time 2h	-4.82222222666667	4.61941688316091	-1.04390280172483	0.32704543157344
PARP1KO vs. ALC1KO in Time 2h	6.54444444666666	4.61941688316092	1.41672523008759	0.194305008686909
PARP1KO vs. WT in Time 2h	-13.38888889	4.61941688316092	-2.89839372125219	0.0199407600005575
PARP1KO vs. ALC1KOPARP1KO in Time 2h	1.72222222	4.61941688316092	0.372822428362763	0.718964761240989
ALC1KOPARP1KO vs. PARP1KO in Time 2h	-1.72222222	4.61941688316091	-0.372822428362763	0.71896476124099
ALC1KOPARP1KO vs. ALC1KO in Time 2h	4.82222222666667	4.61941688316092	1.04390280172483	0.32704543157344
ALC1KOPARP1KO vs. WT in Time 2h	-15.11111111	4.61941688316091	-3.27121614961496	0.0113355916961444
	Estimate	Std. Error	t value	Pr(>|t|)
WT vs. ALC1KO in Time 4h	25.6555555566666	2.94700098560667	8.70564878735025	2.36474253843872e-05
WT vs. PARP1KO in Time 4h	14.1333333366667	2.94700098560667	4.79583597212716	0.00136277121174079
WT vs. ALC1KOPARP1KO in Time 4h	12.4666666666667	2.94700098560667	4.23028927630314	0.00287522676347967
ALC1KO vs. WT in Time 4h	-25.6555555566667	2.94700098560666	-8.70564878735025	2.3647425384387e-05
ALC1KO vs. PARP1KO in Time 4h	-11.52222222	2.94700098560666	-3.90981281522309	0.0044822300200186
ALC1KO vs. ALC1KOPARP1KO in Time 4h	-13.18888889	2.94700098560666	-4.47535951104711	0.00206847786110044
PARP1KO vs. ALC1KO in Time 4h	11.52222222	2.94700098560667	3.90981281522308	0.00448223002001862
PARP1KO vs. WT in Time 4h	-14.1333333366667	2.94700098560667	-4.79583597212717	0.00136277121174079
PARP1KO vs. ALC1KOPARP1KO in Time 4h	-1.66666667	2.94700098560667	-0.565546695824028	0.587208761244944
ALC1KOPARP1KO vs. PARP1KO in Time 4h	1.66666667	2.94700098560666	0.565546695824028	0.587208761244943
ALC1KOPARP1KO vs. ALC1KO in Time 4h	13.18888889	2.94700098560667	4.47535951104711	0.00206847786110044
ALC1KOPARP1KO vs. WT in Time 4h	-12.4666666666667	2.94700098560666	-4.23028927630314	0.00287522676347967
	Estimate	Std. Error	t value	Pr(>|t|)
WT vs. ALC1KO in Time 6h	34.61111111	3.55089103646649	9.747162262811	1.02732215316808e-05
WT vs. PARP1KO in Time 6h	19.5222222233333	3.55089103646649	5.49783759142325	0.000575210034763358
WT vs. ALC1KOPARP1KO in Time 6h	19.9666666666667	3.55089103646649	5.62300179352606	0.000496757913056301
ALC1KO vs. WT in Time 6h	-34.61111111	3.55089103646648	-9.74716226281101	1.02732215316807e-05
ALC1KO vs. PARP1KO in Time 6h	-15.0888888866667	3.55089103646648	-4.24932467138775	0.00280172251246753
ALC1KO vs. ALC1KOPARP1KO in Time 6h	-14.6444444433333	3.55089103646648	-4.12416046928495	0.00332508800054009
PARP1KO vs. ALC1KO in Time 6h	15.0888888866667	3.55089103646649	4.24932467138775	0.00280172251246755
PARP1KO vs. WT in Time 6h	-19.5222222233333	3.55089103646649	-5.49783759142326	0.000575210034763353
PARP1KO vs. ALC1KOPARP1KO in Time 6h	0.44444444333333	3.55089103646649	0.125164202102805	0.903481679641302
ALC1KOPARP1KO vs. PARP1KO in Time 6h	-0.444444443333329	3.55089103646648	-0.125164202102805	0.903481679641302
ALC1KOPARP1KO vs. ALC1KO in Time 6h	14.6444444433333	3.55089103646649	4.12416046928495	0.0033250880005401
ALC1KOPARP1KO vs. WT in Time 6h	-19.9666666666667	3.55089103646648	-5.62300179352606	0.000496757913056298

write.table(output, file = "Figure3D_Stats_New_All.txt", quote = F, sep = "\t", row.names = T, col.names = NA)